program mksurfdat

!-----------------------------------------------------------------------
!BOP
!
! !PROGRAM: mksurfdat
!
! !DESCRIPTION:
! Creates land model surface dataset from original "raw" data files.
! Surface dataset contains model grid, pfts, inland water, glacier,
! soil texture, soil color, LAI and SAI, urban fraction, and urban
! parameters.
!
! !USES:
    use shr_kind_mod       , only : r8 => shr_kind_r8, r4 => shr_kind_r4
    use fileutils          , only : opnfil, getavu
    use mklaiMod           , only : mklai
    use mkpctPftTypeMod    , only : pct_pft_type, get_pct_p2l_array, get_pct_l2g_array, update_max_array
    use mkpftConstantsMod  , only : natpft_lb, natpft_ub, cft_lb, cft_ub, num_cft
    use mkpftMod           , only : pft_idx, pft_frc, mkpft, mkpftInit, mkpft_parse_oride
    use mksoilMod          , only : soil_sand, soil_clay, mksoiltex, mksoilInit, &
                                    soil_color, mksoilcol, mkorganic, &
                                    soil_fmax, mkfmax
    use mkvocefMod         , only : mkvocef
    use mklanwatMod        , only : mklakwat, mkwetlnd, mklakparams
    use mkglacierregionMod , only : mkglacierregion
    use mkglcmecMod        , only : nglcec, mkglcmec, mkglcmecInit, mkglacier
    use mkharvestMod       , only : mkharvest, mkharvest_init, mkharvest_fieldname
    use mkharvestMod       , only : mkharvest_numtypes, mkharvest_parse_oride
    use mkharvestMod       , only : harvestDataType
    use mkurbanparCommonMod, only : mkelev
    use mkurbanparMod      , only : mkurbanInit, mkurban, mkurbanpar, numurbl
    use mkutilsMod         , only : normalize_classes_by_gcell
    use mkfileMod          , only : mkfile
    use mkvarpar           , only : nlevsoi, elev_thresh, numstdpft
    use mkvarctl
    use nanMod             , only : nan, bigint
    use mkncdio            , only : check_ret, ncd_put_time_slice
    use mkdomainMod        , only : domain_type, domain_read_map, domain_read, &
                                    domain_write
    use mkgdpMod           , only : mkgdp
    use mkpeatMod          , only : mkpeat
    use mksoildepthMod          , only : mksoildepth
    use mkagfirepkmonthMod , only : mkagfirepkmon
    use mktopostatsMod     , only : mktopostats
    use mkVICparamsMod     , only : mkVICparams
!
! !ARGUMENTS:
    implicit none

    include 'netcdf.inc'
!
! !REVISION HISTORY:
! Authors: Gordon Bonan, Sam Levis and Mariana Vertenstein
! Revised: Nan Rosenbloom to add fmax processing.
! 3/18/08: David Lawrence added organic matter processing
! 1/22/09: Keith Oleson added urban parameter processing
! 2/11/13: Sam Levis added abm, peat, and gdp processing for new fire model
!
!
! !LOCAL VARIABLES:
!EOP
    integer  :: nsoicol                     ! number of model color classes
    integer  :: k,m,n                       ! indices
    integer  :: ni,nj,ns_o                  ! indices
    integer  :: ier                         ! error status
    integer  :: ndiag,nfdyn                 ! unit numbers
    integer  :: ncid                        ! netCDF id
    integer  :: omode                       ! netCDF output mode
    integer  :: varid                       ! netCDF variable id
    integer  :: ret                         ! netCDF return status
    integer  :: ntim                        ! time sample for dynamic land use
    integer  :: year                        ! year for dynamic land use
    integer  :: year2                       ! year for dynamic land use for harvest file
    logical  :: all_veg                     ! if gridcell will be 100% vegetated land-cover
    real(r8) :: suma                        ! sum for error check
    character(len=256) :: fgrddat           ! grid data file
    character(len=256) :: fsurdat           ! output surface data file name (if blank, do not output a surface dataset)
    character(len=256) :: fsurlog           ! output surface log file name
    character(len=256) :: fdyndat           ! dynamic landuse data file name
    character(len=256) :: fname             ! generic filename
    character(len=256) :: fhrvname          ! generic harvest filename
    character(len=256) :: string            ! string read in
    integer  :: t1                          ! timer
    real(r8),parameter :: p5  = 0.5_r8      ! constant
    real(r8),parameter :: p25 = 0.25_r8     ! constant

    real(r8), allocatable  :: landfrac_pft(:)    ! PFT data: % land per gridcell
    real(r8), allocatable  :: pctlnd_pft(:)      ! PFT data: % of gridcell for PFTs
    real(r8), allocatable  :: pctlnd_pft_dyn(:)  ! PFT data: % of gridcell for dyn landuse PFTs
    integer , allocatable  :: pftdata_mask(:)    ! mask indicating real or fake land type
    type(pct_pft_type), allocatable :: pctnatpft(:)     ! % of grid cell that is nat veg, and breakdown into PFTs
    type(pct_pft_type), allocatable :: pctnatpft_max(:) ! % of grid cell maximum PFTs of the time series
    type(pct_pft_type), allocatable :: pctcft(:)        ! % of grid cell that is crop, and breakdown into CFTs
    type(pct_pft_type), allocatable :: pctcft_max(:)    ! % of grid cell maximum CFTs of the time series
    real(r8)               :: harvest_initval    ! initial value for harvest variables
    real(r8), pointer      :: harvest1D(:)       ! harvest 1D data: normalized harvesting
    real(r8), pointer      :: harvest2D(:,:)     ! harvest 1D data: normalized harvesting
    real(r8), allocatable  :: pctgla(:)          ! percent of grid cell that is glacier  
    real(r8), allocatable  :: pctglc_gic(:)      ! percent of grid cell that is gic (% of glc landunit)
    real(r8), allocatable  :: pctglc_icesheet(:) ! percent of grid cell that is ice sheet (% of glc landunit)
    real(r8), allocatable  :: pctglcmec(:,:)     ! glacier_mec pct coverage in each class (% of landunit)
    real(r8), allocatable  :: topoglcmec(:,:)    ! glacier_mec sfc elevation in each gridcell and class
    real(r8), allocatable  :: pctglcmec_gic(:,:) ! GIC pct coverage in each class (% of landunit)
    real(r8), allocatable  :: pctglcmec_icesheet(:,:) ! icesheet pct coverage in each class (% of landunit)
    real(r8), allocatable  :: elevclass(:)       ! glacier_mec elevation classes
    integer,  allocatable  :: glacier_region(:)  ! glacier region ID
    real(r8), allocatable  :: pctlak(:)          ! percent of grid cell that is lake
    real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
    real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
    real(r8), allocatable  :: urbn_classes(:,:)  ! percent cover of each urban class, as % of total urban area
    real(r8), allocatable  :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell
    real(r8), allocatable  :: elev(:)            ! glc elevation (m)
    real(r8), allocatable  :: fmax(:)            ! fractional saturated area
    integer , allocatable  :: soicol(:)          ! soil color                            
    real(r8), allocatable  :: pctsand(:,:)       ! soil texture: percent sand            
    real(r8), allocatable  :: pctclay(:,:)       ! soil texture: percent clay            
    real(r8), allocatable  :: ef1_btr(:)         ! Isoprene emission factor for broadleaf
    real(r8), allocatable  :: ef1_fet(:)         ! Isoprene emission factor for fine/everg
    real(r8), allocatable  :: ef1_fdt(:)         ! Isoprene emission factor for fine/dec
    real(r8), allocatable  :: ef1_shr(:)         ! Isoprene emission factor for shrubs
    real(r8), allocatable  :: ef1_grs(:)         ! Isoprene emission factor for grasses
    real(r8), allocatable  :: ef1_crp(:)         ! Isoprene emission factor for crops
    real(r8), allocatable  :: organic(:,:)       ! organic matter density (kg/m3)            
    real(r8), allocatable  :: gdp(:)             ! GDP (x1000 1995 US$/capita)
    real(r8), allocatable  :: fpeat(:)           ! peatland fraction of gridcell
    real(r8), allocatable  :: soildepth(:)       ! soil depth (m)
    integer , allocatable  :: agfirepkmon(:)     ! agricultural fire peak month
    integer , allocatable  :: urban_region(:)    ! urban region ID
    real(r8), allocatable  :: topo_stddev(:)     ! standard deviation of elevation (m)
    real(r8), allocatable  :: slope(:)           ! topographic slope (degrees)
    real(r8), allocatable  :: vic_binfl(:)       ! VIC b parameter (unitless)
    real(r8), allocatable  :: vic_ws(:)          ! VIC Ws parameter (unitless)
    real(r8), allocatable  :: vic_dsmax(:)       ! VIC Dsmax parameter (mm/day)
    real(r8), allocatable  :: vic_ds(:)          ! VIC Ds parameter (unitless)
    real(r8), allocatable  :: lakedepth(:)       ! lake depth (m)

    real(r8) :: std_elev = -999.99_r8            ! Standard deviation of elevation (m) to use for entire grid

    integer, allocatable :: harvind1D(:)         ! Indices of 1D harvest fields
    integer, allocatable :: harvind2D(:)         ! Indices of 2D harvest fields

    ! NOTE(bja, 2015-01) added to work around a ?bug? causing 1x1_urbanc_alpha to abort. See
    !/glade/p/cesm/cseg/inputdata/lnd/clm2/surfdata_map/README_c141219
    logical :: urban_skip_abort_on_invalid_data_check

    type(domain_type) :: ldomain

    character(len=32) :: subname = 'mksrfdat'  ! program name
    type(harvestDataType) :: harvdata

    namelist /clmexp/              &
	 mksrf_fgrid,              &	
	 mksrf_gridtype,           &	
         mksrf_fvegtyp,            &
         mksrf_fhrvtyp,            &
	 mksrf_fsoitex,            &
         mksrf_forganic,           &
         mksrf_fsoicol,            &
         mksrf_fvocef,             &
         mksrf_flakwat,            &
         mksrf_fwetlnd,            &
         mksrf_fglacier,           &
         mksrf_fglacierregion,     &
         mksrf_furbtopo,           &
         mksrf_fmax,               &
         mksrf_furban,             &
         mksrf_flai,               &
         mksrf_fdynuse,            &
         mksrf_fgdp,               &
         mksrf_fpeat,              &
         mksrf_fsoildepth,         &
         mksrf_fabm,               &
         mksrf_ftopostats,         &
         mksrf_fvic,               &
         mksrf_fch4,               &
         nglcec,                   &
         numpft,                   &
         soil_color,               &
         soil_sand,                &
         soil_fmax,                &
         soil_clay,                &
         pft_idx,                  &
         all_veg,                  &
         pft_frc,                  &
         all_urban,                &
         no_inlandwet,             &
         map_fpft,                 &
         map_flakwat,              &
         map_fwetlnd,              &
         map_fglacier,             &
         map_fglacierregion,       &
         map_fsoitex,              &
         map_fsoicol,              &
         map_furban,               &
         map_furbtopo,             &
         map_fmax,                 &
         map_forganic,             &
         map_fvocef,               &
         map_flai,                 &
         map_fharvest,             &
         map_fgdp,                 &
         map_fpeat,                &
         map_fsoildepth,           &
         map_fabm,                 &
         map_ftopostats,           &
         map_fvic,                 &
         map_fch4,                 &
         gitdescribe,              &
         outnc_large_files,        &
         outnc_double,             &
         outnc_dims,               &
         outnc_vic,                &
         outnc_3dglc,              &
         fsurdat,                  &
         fdyndat,                  &   
         fsurlog,                  &
         std_elev,                 &
         urban_skip_abort_on_invalid_data_check

!-----------------------------------------------------------------------

    ! ======================================================================
    ! Read input namelist
    ! ======================================
    ! Must specify settings for the output grid:
    ! ======================================
    !    mksrf_fgrid -- Grid dataset
    ! ======================================
    ! Must specify settings for input high resolution datafiles
    ! ======================================
    !    mksrf_fglacier - Glacier dataset
    !    mksrf_fglacierregion - Glacier region ID dataset
    !    mksrf_flai ----- Leaf Area Index dataset
    !    mksrf_flakwat -- Lake water dataset
    !    mksrf_fwetlnd -- Wetland water dataset
    !    mksrf_forganic - Organic soil carbon dataset
    !    mksrf_fmax ----- Max fractional saturated area dataset
    !    mksrf_fsoicol -- Soil color dataset
    !    mksrf_fsoitex -- Soil texture dataset
    !    mksrf_furbtopo-- Topography dataset (for limiting urban areas)
    !    mksrf_furban --- Urban dataset
    !    mksrf_fvegtyp -- PFT vegetation type dataset
    !    mksrf_fhrvtyp -- harvest type dataset
    !    mksrf_fvocef  -- Volatile Organic Compund Emission Factor dataset
    !    mksrf_fgdp ----- GDP dataset
    !    mksrf_fpeat ---- Peatland dataset
    !    mksrf_fsoildepth Soil depth dataset
    !    mksrf_fabm ----- Agricultural fire peak month dataset
    !    mksrf_ftopostats Topography statistics dataset
    !    mksrf_fvic ----- VIC parameters dataset
    !    mksrf_fch4 ----- inversion-derived CH4 parameters dataset
    ! ======================================
    ! Must specify mapping file for the different datafiles above
    ! ======================================
    !    map_fpft -------- Mapping for mksrf_fvegtyp
    !    map_flakwat ----- Mapping for mksrf_flakwat
    !    map_fwetlnd ----- Mapping for mksrf_fwetlnd
    !    map_fglacier ---- Mapping for mksrf_fglacier
    !    map_fglacierregion - Mapping for mksrf_fglacierregion
    !    map_fsoitex ----- Mapping for mksrf_fsoitex
    !    map_fsoicol ----- Mapping for mksrf_fsoicol
    !    map_furban ------ Mapping for mksrf_furban
    !    map_furbtopo ---- Mapping for mksrf_furbtopo
    !    map_fmax -------- Mapping for mksrf_fmax
    !    map_forganic ---- Mapping for mksrf_forganic
    !    map_fvocef ------ Mapping for mksrf_fvocef
    !    map_flai -------- Mapping for mksrf_flai
    !    map_fharvest ---- Mapping for mksrf_flai harvesting
    !    map_fgdp -------- Mapping for mksrf_fgdp
    !    map_fpeat ------- Mapping for mksrf_fpeat
    !    map_fsoildepth -- Mapping for mksrf_fsoildepth
    !    map_fabm -------- Mapping for mksrf_fabm
    !    map_ftopostats -- Mapping for mksrf_ftopostats
    !    map_fvic -------- Mapping for mksrf_fvic
    !    map_fch4 -------- Mapping for mksrf_fch4
    ! ======================================
    ! Optionally specify setting for:
    ! ======================================
    !    mksrf_fdynuse ----- ASCII text file that lists each year of pft files to use
    !    mksrf_gridtype ---- Type of grid (default is 'global')
    !    outnc_double ------ If output should be in double precision
    !    outnc_large_files - If output should be in NetCDF large file format
    !    outnc_vic --------- Output fields needed for VIC
    !    outnc_3dglc ------- Output 3D glacier fields (normally only needed for comparasion)
    !    nglcec ------------ If you want to change the number of Glacier elevation classes
    !    gitdescribe ------- Description of this version from git
    ! ======================================
    ! Optional settings to change values for entire area
    ! ======================================
    !    all_urban --------- If entire area is urban
    !    all_veg ----------- If entire area is to be vegetated (pft_idx and pft_frc then required)
    !    no_inlandwet ------ If wetland should be set to 0% over land
    !    soil_color -------- If you want to change the soil_color to this value everywhere
    !    soil_clay --------- If you want to change the soil_clay % to this value everywhere
    !    soil_fmax --------- If you want to change the soil_fmax  to this value everywhere
    !    soil_sand --------- If you want to change the soil_sand % to this value everywhere
    !    pft_idx ----------- If you want to change to 100% veg covered with given PFT indices
    !    pft_frc ----------- Fractions that correspond to the pft_idx above
    ! ==================
    !    numpft            (if different than default of 16)
    ! ======================================
    ! Optional settings to work around urban bug?
    ! ======================================
    !    urban_skip_abort_on_invalid_data_check
    ! ======================================================================

    write(6,*) 'Attempting to initialize control settings .....'

    mksrf_gridtype    = 'global'
    outnc_large_files = .false.
    outnc_double      = .true.
    outnc_vic         = .false.
    outnc_3dglc       = .false.
    all_urban         = .false.
    all_veg           = .false.
    no_inlandwet      = .true.

    ! default value for bug work around
    urban_skip_abort_on_invalid_data_check = .false.

    read(5, clmexp, iostat=ier)
    if (ier /= 0) then
       write(6,*)'error: namelist input resulted in error code ',ier
       call abort()
    endif

    write (6,*) 'Attempting to create surface boundary data .....'
    write (6,'(72a1)') ("-",n=1,60)

    ! ----------------------------------------------------------------------
    ! Error check namelist input
    ! ----------------------------------------------------------------------
    
    if (urban_skip_abort_on_invalid_data_check) then
       write(6, *) "WARNING: aborting on invalid data check in urban has been disabled!"
       write(6, *) "WARNING: urban data may be invalid!"
    end if
    
    if (mksrf_fgrid /= ' ')then
       fgrddat = mksrf_fgrid
       write(6,*)'mksrf_fgrid = ',mksrf_fgrid
    else
       write (6,*)'must specify mksrf_fgrid'
       call abort()
    endif

    if (trim(mksrf_gridtype) == 'global' .or. &
        trim(mksrf_gridtype) == 'regional') then
       write(6,*)'mksrf_gridtype = ',trim(mksrf_gridtype)
    else
       write(6,*)'mksrf_gridtype = ',trim(mksrf_gridtype)
       write (6,*)'illegal mksrf_gridtype, must be global or regional '
       call abort()
    endif
    if ( outnc_large_files )then
       write(6,*)'Output file in NetCDF 64-bit large_files format'
    end if
    if ( outnc_double )then
       write(6,*)'Output ALL data in file as 64-bit'
    end if
    if ( outnc_vic )then
       write(6,*)'Output VIC fields'
    end if
    if ( outnc_3dglc )then
       write(6,*)'Output optional 3D glacier fields (mostly used for verification of the glacier model)'
    end if
    if ( outnc_3dglc )then
       write(6,*)'Output optional 3D glacier fields (mostly used for verification of the glacier model)'
    end if
    if ( all_urban )then
       write(6,*) 'Output ALL data in file as 100% urban'
    end if
    if ( no_inlandwet )then
       write(6,*) 'Set wetland to 0% over land'
    end if
    if (nglcec <= 0) then
       write(6,*) 'nglcec must be at least 1'
       call abort()
    end if

    !
    ! Call module initialization routines
    !
    call mksoilInit( )
    call mkpftInit( zero_out_l=all_urban, all_veg_l=all_veg )
    allocate ( elevclass(nglcec+1) )
    call mkglcmecInit (elevclass)
    call mkurbanInit (mksrf_furban)

    if ( all_veg )then
       write(6,*) 'Output ALL data in file as 100% vegetated'
    end if

    ! ----------------------------------------------------------------------
    ! Determine land model grid, fractional land and land mask
    ! ----------------------------------------------------------------------
    
    write(6,*)'calling domain_read'
    if ( .not. domain_read_map(ldomain, fgrddat) )then
        call domain_read(ldomain, fgrddat)
    end if
    write(6,*)'finished domain_read'
    
    ! Invalidate mask and frac for ldomain 

    !ldomain%mask = bigint
    !ldomain%frac = nan

    ! Determine if will have 1d output

    if (ldomain%ni /= -9999 .and. ldomain%nj /= -9999) then
       write(6,*)'fsurdat is 2d lat/lon grid'
       write(6,*)'nlon= ',ldomain%ni,' nlat= ',ldomain%nj
       if (outnc_dims == 1) then
          write(6,*)' writing output file in 1d gridcell format'
       end if
    else
       write(6,*)'fsurdat is 1d gridcell grid'
       outnc_dims = 1
    end if

    outnc_1d = .false.
    if ((ldomain%ni == -9999 .and. ldomain%nj == -9999) .or. outnc_dims==1) then
       outnc_1d = .true.
       write(6,*)'output file will be 1d'
    end if

    ! ----------------------------------------------------------------------
    ! Allocate and initialize dynamic memory
    ! ----------------------------------------------------------------------

    ns_o = ldomain%ns
    allocate ( landfrac_pft(ns_o)                 , &
               pctlnd_pft(ns_o)                   , & 
               pftdata_mask(ns_o)                 , & 
               pctnatpft(ns_o)                    , &
               pctnatpft_max(ns_o)                , &
               pctcft(ns_o)                       , &
               pctcft_max(ns_o)                   , &
               pctgla(ns_o)                       , & 
               pctlak(ns_o)                       , & 
               pctwet(ns_o)                       , & 
               pcturb(ns_o)                       , &
               urban_region(ns_o)                 , &
               urbn_classes(ns_o,numurbl)         , &
               urbn_classes_g(ns_o,numurbl)       , &
               pctsand(ns_o,nlevsoi)              , & 
               pctclay(ns_o,nlevsoi)              , & 
               soicol(ns_o)                       , & 
               gdp(ns_o)                          , & 
               fpeat(ns_o)                        , & 
               soildepth(ns_o)                    , & 
               agfirepkmon(ns_o)                  , & 
               topo_stddev(ns_o)                  , &
               slope(ns_o)                        , &
               vic_binfl(ns_o)                    , &
               vic_ws(ns_o)                       , &
               vic_dsmax(ns_o)                    , &
               vic_ds(ns_o)                       , &
               lakedepth(ns_o)                    , &
               glacier_region(ns_o)               )       
    landfrac_pft(:)       = spval 
    pctlnd_pft(:)         = spval
    pftdata_mask(:)       = -999
    pctgla(:)             = spval
    pctlak(:)             = spval
    pctwet(:)             = spval
    pcturb(:)             = spval
    urban_region(:)       = -999
    urbn_classes(:,:)     = spval
    urbn_classes_g(:,:)   = spval
    pctsand(:,:)          = spval
    pctclay(:,:)          = spval
    soicol(:)             = -999
    gdp(:)                = spval
    fpeat(:)              = spval
    soildepth(:)          = spval
    agfirepkmon(:)        = -999
    topo_stddev(:)        = spval
    slope(:)              = spval
    vic_binfl(:)          = spval
    vic_ws(:)             = spval
    vic_dsmax(:)          = spval
    vic_ds(:)             = spval
    lakedepth(:)          = spval
    glacier_region(:)     = -999

    ! ----------------------------------------------------------------------
    ! Open diagnostic output log file
    ! ----------------------------------------------------------------------
    
    if (fsurlog == ' ') then
       write(6,*)' must specify fsurlog in namelist'
       stop
    else
       ndiag = getavu(); call opnfil (fsurlog, ndiag, 'f')
    end if
    
    if (urban_skip_abort_on_invalid_data_check) then
       write(ndiag, *) "WARNING: aborting on invalid data check in urban has been disabled!"
       write(ndiag, *) "WARNING: urban data may be invalid!"
    end if
    
    if (mksrf_fgrid /= ' ')then
       write (ndiag,*)'using fractional land data from file= ', &
            trim(mksrf_fgrid),' to create the surface dataset'
    endif

    if (trim(mksrf_gridtype) == 'global' .or. &
        trim(mksrf_gridtype) == 'regional') then
       write(6,*)'mksrf_gridtype = ',trim(mksrf_gridtype)
    endif

    write(ndiag,*) 'PFTs from:                   ',trim(mksrf_fvegtyp)
    write(ndiag,*) 'harvest from:                ',trim(mksrf_fhrvtyp)
    write(ndiag,*) 'fmax from:                   ',trim(mksrf_fmax)
    write(ndiag,*) 'glaciers from:               ',trim(mksrf_fglacier)
    write(ndiag,*) '           with:             ', nglcec, ' glacier elevation classes'
    write(ndiag,*) 'glacier region ID from:      ',trim(mksrf_fglacierregion)
    write(ndiag,*) 'urban topography from:       ',trim(mksrf_furbtopo)
    write(ndiag,*) 'urban from:                  ',trim(mksrf_furban)
    write(ndiag,*) 'inland lake from:            ',trim(mksrf_flakwat)
    write(ndiag,*) 'inland wetland from:         ',trim(mksrf_fwetlnd)
    write(ndiag,*) 'soil texture from:           ',trim(mksrf_fsoitex)
    write(ndiag,*) 'soil organic from:           ',trim(mksrf_forganic)
    write(ndiag,*) 'soil color from:             ',trim(mksrf_fsoicol)
    write(ndiag,*) 'VOC emission factors from:   ',trim(mksrf_fvocef)
    write(ndiag,*) 'gdp from:                    ',trim(mksrf_fgdp)
    write(ndiag,*) 'peat from:                   ',trim(mksrf_fpeat)
    write(ndiag,*) 'soil depth from:             ',trim(mksrf_fsoildepth)
    write(ndiag,*) 'abm from:                    ',trim(mksrf_fabm)
    write(ndiag,*) 'topography statistics from:  ',trim(mksrf_ftopostats)
    write(ndiag,*) 'VIC parameters from:         ',trim(mksrf_fvic)
    write(ndiag,*) 'CH4 parameters from:         ',trim(mksrf_fch4)
    write(ndiag,*)' mapping for pft              ',trim(map_fpft)
    write(ndiag,*)' mapping for lake water       ',trim(map_flakwat)
    write(ndiag,*)' mapping for wetland          ',trim(map_fwetlnd)
    write(ndiag,*)' mapping for glacier          ',trim(map_fglacier)
    write(ndiag,*)' mapping for glacier region   ',trim(map_fglacierregion)
    write(ndiag,*)' mapping for soil texture     ',trim(map_fsoitex)
    write(ndiag,*)' mapping for soil color       ',trim(map_fsoicol)
    write(ndiag,*)' mapping for soil organic     ',trim(map_forganic)
    write(ndiag,*)' mapping for urban            ',trim(map_furban)
    write(ndiag,*)' mapping for fmax             ',trim(map_fmax)
    write(ndiag,*)' mapping for VOC pct emis     ',trim(map_fvocef)
    write(ndiag,*)' mapping for harvest          ',trim(map_fharvest)
    write(ndiag,*)' mapping for lai/sai          ',trim(map_flai)
    write(ndiag,*)' mapping for urb topography   ',trim(map_furbtopo)
    write(ndiag,*)' mapping for GDP              ',trim(map_fgdp)
    write(ndiag,*)' mapping for peatlands        ',trim(map_fpeat)
    write(ndiag,*)' mapping for soil depth       ',trim(map_fsoildepth)
    write(ndiag,*)' mapping for ag fire pk month ',trim(map_fabm)
    write(ndiag,*)' mapping for topography stats ',trim(map_ftopostats)
    write(ndiag,*)' mapping for VIC parameters   ',trim(map_fvic)
    write(ndiag,*)' mapping for CH4 parameters   ',trim(map_fch4)

    if (mksrf_fdynuse /= ' ') then
       write(6,*)'mksrf_fdynuse = ',trim(mksrf_fdynuse)
    end if

    ! ----------------------------------------------------------------------
    ! Make surface dataset fields
    ! ----------------------------------------------------------------------

    ! Make PFTs [pctnatpft, pctcft] from dataset [fvegtyp]

    call mkpft(ldomain, mapfname=map_fpft, fpft=mksrf_fvegtyp, &
         ndiag=ndiag, pctlnd_o=pctlnd_pft, pctnatpft_o=pctnatpft, pctcft_o=pctcft)

    ! Create harvesting data at model resolution
    if (all_veg) then
       ! In this case, we don't call mkharvest, so we want the harvest variables to be
       ! initialized reasonably.
       harvest_initval = 0._r8
    else
       harvest_initval = spval
    end if
    call mkharvest_init( ns_o, harvest_initval, harvdata, mksrf_fhrvtyp )
    if ( .not. all_veg )then

       call mkharvest( ldomain, mapfname=map_fharvest, datfname=mksrf_fhrvtyp, &
                       ndiag=ndiag, harvdata=harvdata )
    end if

    ! Make inland water [pctlak, pctwet] [flakwat] [fwetlnd]

    call mklakwat (ldomain, mapfname=map_flakwat, datfname=mksrf_flakwat, &
         ndiag=ndiag, zero_out=all_urban.or.all_veg, lake_o=pctlak)

    call mkwetlnd (ldomain, mapfname=map_fwetlnd, datfname=mksrf_fwetlnd, &
         ndiag=ndiag, zero_out=all_urban.or.all_veg.or.no_inlandwet, swmp_o=pctwet)

    ! Make glacier fraction [pctgla] from [fglacier] dataset

    call mkglacier (ldomain, mapfname=map_fglacier, datfname=mksrf_fglacier, &
         ndiag=ndiag, zero_out=all_urban.or.all_veg, glac_o=pctgla)

    ! Make glacier region ID [glacier_region] from [fglacierregion] dataset

    call mkglacierregion (ldomain, mapfname=map_fglacierregion, &
         datfname=mksrf_fglacierregion, ndiag=ndiag, &
         glacier_region_o = glacier_region)

    ! Make soil texture [pctsand, pctclay]  [fsoitex]

    call mksoiltex (ldomain, mapfname=map_fsoitex, datfname=mksrf_fsoitex, &
         ndiag=ndiag, sand_o=pctsand, clay_o=pctclay)
    ! Make soil color classes [soicol] [fsoicol]

    call mksoilcol (ldomain, mapfname=map_fsoicol, datfname=mksrf_fsoicol, &
         ndiag=ndiag, soil_color_o=soicol, nsoicol=nsoicol)

    ! Make fmax [fmax] from [fmax] dataset

    allocate(fmax(ns_o))
    fmax(:) = spval
    call mkfmax (ldomain, mapfname=map_fmax, datfname=mksrf_fmax, &
         ndiag=ndiag, fmax_o=fmax)

    ! Make GDP data [gdp] from [gdp]

    call mkgdp (ldomain, mapfname=map_fgdp, datfname=mksrf_fgdp, &
         ndiag=ndiag, gdp_o=gdp)

    ! Make peat data [fpeat] from [peatf]

    call mkpeat (ldomain, mapfname=map_fpeat, datfname=mksrf_fpeat, &
         ndiag=ndiag, peat_o=fpeat)

    ! Make soil depth data [soildepth] from [soildepthf]

    call mksoildepth (ldomain, mapfname=map_fsoildepth, datfname=mksrf_fsoildepth, &
         ndiag=ndiag, soildepth_o=soildepth)

    ! Make agricultural fire peak month data [abm] from [abm]

    call mkagfirepkmon (ldomain, mapfname=map_fabm, datfname=mksrf_fabm, &
         ndiag=ndiag, agfirepkmon_o=agfirepkmon)

    ! Make urban fraction [pcturb] from [furban] dataset

    call mkurban (ldomain, mapfname=map_furban, datfname=mksrf_furban, &
         ndiag=ndiag, zero_out=all_veg, urbn_o=pcturb, urbn_classes_o=urbn_classes, &
         region_o=urban_region)

    ! Make elevation [elev] from [ftopo, ffrac] dataset
    ! Used only to screen pcturb  
    ! Screen pcturb by elevation threshold from elev dataset

    if ( .not. all_urban .and. .not. all_veg )then
       allocate(elev(ns_o))
       elev(:) = spval
       ! NOTE(wjs, 2016-01-15) This uses the 'TOPO_ICE' variable for historical reasons
       ! (this same dataset used to be used for glacier-related purposes as well).
       ! TODO(wjs, 2016-01-15) A better solution for this urban screening would probably
       ! be to modify the raw urban data; in that case, I believe we could remove
       ! furbtopo.
       call mkelev (ldomain, mapfname=map_furbtopo, datfname=mksrf_furbtopo, &
         varname='TOPO_ICE', ndiag=ndiag, elev_o=elev)

       where (elev .gt. elev_thresh)
         pcturb = 0._r8
       end where
       deallocate(elev)
    end if
    
    ! Compute topography statistics [topo_stddev, slope] from [ftopostats]
    call mktopostats (ldomain, mapfname=map_ftopostats, datfname=mksrf_ftopostats, &
         ndiag=ndiag, topo_stddev_o=topo_stddev, slope_o=slope, std_elev=std_elev)

    ! Make VIC parameters [binfl, ws, dsmax, ds] from [fvic]
    if ( outnc_vic )then
       call mkVICparams (ldomain, mapfname=map_fvic, datfname=mksrf_fvic, ndiag=ndiag, &
            binfl_o=vic_binfl, ws_o=vic_ws, dsmax_o=vic_dsmax, ds_o=vic_ds)
    end if

    ! Make lake depth [lakedepth] from [flakwat]
    call mklakparams (ldomain, mapfname=map_flakwat, datfname=mksrf_flakwat, ndiag=ndiag, &
         lakedepth_o=lakedepth)

    ! Make organic matter density [organic] [forganic]
    allocate (organic(ns_o,nlevsoi))
    organic(:,:) = spval
    call mkorganic (ldomain, mapfname=map_forganic, datfname=mksrf_forganic, &
         ndiag=ndiag, organic_o=organic)

    ! Make VOC emission factors for isoprene &
    ! [ef1_btr,ef1_fet,ef1_fdt,ef1_shr,ef1_grs,ef1_crp]

    allocate ( ef1_btr(ns_o) , & 
               ef1_fet(ns_o) , & 
               ef1_fdt(ns_o) , & 
               ef1_shr(ns_o) , & 
               ef1_grs(ns_o) , & 
               ef1_crp(ns_o) )
    ef1_btr(:) = 0._r8
    ef1_fet(:) = 0._r8
    ef1_fdt(:) = 0._r8
    ef1_shr(:) = 0._r8
    ef1_grs(:) = 0._r8
    ef1_crp(:) = 0._r8
    
    call mkvocef (ldomain, mapfname=map_fvocef, datfname=mksrf_fvocef, ndiag=ndiag, &
         ef_btr_o=ef1_btr, ef_fet_o=ef1_fet, ef_fdt_o=ef1_fdt,  &
         ef_shr_o=ef1_shr, ef_grs_o=ef1_grs, ef_crp_o=ef1_crp)

    ! Do landuse changes such as for the poles, etc.

    call change_landuse( ldomain, dynpft=.false. )

    do n = 1,ns_o

       ! Assume wetland, glacier and/or lake when dataset landmask implies ocean 
       ! (assume medium soil color (15) and loamy texture).
       ! Also set pftdata_mask here

       if (pctlnd_pft(n) < 1.e-6_r8) then
          pftdata_mask(n)  = 0
          soicol(n)        = 15
          if (pctgla(n) < 1.e-6_r8) then
              pctwet(n)    = 100._r8 - pctlak(n)
              pctgla(n)    = 0._r8
          else
              pctwet(n)    = max(100._r8 - pctgla(n) - pctlak(n), 0.0_r8)
          end if
          pcturb(n)        = 0._r8
          call pctnatpft(n)%set_pct_l2g(0._r8)
          call pctcft(n)%set_pct_l2g(0._r8)
          pctsand(n,:)     = 43._r8
          pctclay(n,:)     = 18._r8
          organic(n,:)   = 0._r8
       else
          pftdata_mask(n) = 1
       end if

       ! Truncate all percentage fields on output grid. This is needed to
       ! insure that wt is zero (not a very small number such as
       ! 1e-16) where it really should be zero
       
       do k = 1,nlevsoi
          pctsand(n,k) = float(nint(pctsand(n,k)))
          pctclay(n,k) = float(nint(pctclay(n,k)))
       end do
       pctlak(n) = float(nint(pctlak(n)))
       pctwet(n) = float(nint(pctwet(n)))
       pctgla(n) = float(nint(pctgla(n)))
       
       ! Make sure sum of land cover types does not exceed 100. If it does,
       ! subtract excess from most dominant land cover.
       
       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
       if (suma > 250._r4) then
          write (6,*) subname, ' error: sum of pctlak, pctwet,', &
               'pcturb and pctgla is greater than 250%'
          write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', &
               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n)
          call abort()
       else if (suma > 100._r4) then
          pctlak(n) = pctlak(n) * 100._r8/suma
          pctwet(n) = pctwet(n) * 100._r8/suma
          pcturb(n) = pcturb(n) * 100._r8/suma
          pctgla(n) = pctgla(n) * 100._r8/suma
       end if
       
    end do

    call normalizencheck_landuse(ldomain)

    ! Write out sum of PFT's

    do k = natpft_lb,natpft_ub
       suma = 0._r8
       do n = 1,ns_o
          suma = suma + pctnatpft(n)%get_one_pct_p2g(k)
       enddo
       write(6,*) 'sum over domain of pft ',k,suma
    enddo
    write(6,*)

    do k = cft_lb,cft_ub
       suma = 0._r8
       do n = 1,ns_o
          suma = suma + pctcft(n)%get_one_pct_p2g(k)
       enddo
       write(6,*) 'sum over domain of cft ',k,suma
    enddo
    write(6,*)

    ! Make final values of percent urban by class
    ! This call needs to occur after all corrections are made to pcturb

    call normalize_classes_by_gcell(urbn_classes, pcturb, urbn_classes_g)


    ! Make glacier multiple elevation classes [pctglcmec,topoglcmec] from [fglacier] dataset
    ! This call needs to occur after pctgla has been adjusted for the final time

    allocate (pctglcmec(ns_o,nglcec), &
         topoglcmec(ns_o,nglcec) )
    if ( outnc_3dglc )then
       allocate( &
         pctglcmec_gic(ns_o,nglcec), &
         pctglcmec_icesheet(ns_o,nglcec))
       allocate (pctglc_gic(ns_o))
       allocate (pctglc_icesheet(ns_o))
    end if

    pctglcmec(:,:)          = spval
    topoglcmec(:,:)         = spval

    if ( outnc_3dglc )then
       call mkglcmec (ldomain, mapfname=map_fglacier, &
            datfname_fglacier=mksrf_fglacier, ndiag=ndiag, &
            pctglcmec_o=pctglcmec, topoglcmec_o=topoglcmec, &
            pctglcmec_gic_o=pctglcmec_gic, pctglcmec_icesheet_o=pctglcmec_icesheet, &
            pctglc_gic_o=pctglc_gic, pctglc_icesheet_o=pctglc_icesheet)
    else
       call mkglcmec (ldomain, mapfname=map_fglacier, &
            datfname_fglacier=mksrf_fglacier, ndiag=ndiag, &
            pctglcmec_o=pctglcmec, topoglcmec_o=topoglcmec )
    end if

    ! Determine fractional land from pft dataset

    do n = 1,ns_o
       landfrac_pft(n) = pctlnd_pft(n)/100._r8
    end do

    ! ----------------------------------------------------------------------
    ! Create surface dataset
    ! ----------------------------------------------------------------------

    ! Create netCDF surface dataset.  

    ! If fsurdat is blank, then we do not write a surface dataset - but we may still
    ! write a dynamic landuse file. This is useful if we are creating many datasets at
    ! once, and don't want duplicate surface datasets.
    !
    ! TODO(wjs, 2016-01-26) Ideally, we would also avoid doing the processing of
    ! variables that are just needed by the surface dataset (not by the dynamic landuse
    ! file). However, this would require some analysis of the above code, to determine
    ! which processing is needed (directly or indirectly) in order to create a dynamic
    ! landuse file.

    if (fsurdat /= ' ') then

       call mkfile(ldomain, trim(fsurdat), harvdata, dynlanduse = .false.)

       call domain_write(ldomain, fsurdat)

       call check_ret(nf_open(trim(fsurdat), nf_write, ncid), subname)
       call check_ret(nf_set_fill (ncid, nf_nofill, omode), subname)

       ! Write fields OTHER THAN lai, sai, heights, and urban parameters to netcdf surface dataset

       call check_ret(nf_inq_varid(ncid, 'natpft', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, (/(n,n=natpft_lb,natpft_ub)/)), subname)

       if (num_cft > 0) then
          call check_ret(nf_inq_varid(ncid, 'cft', varid), subname)
          call check_ret(nf_put_var_int(ncid, varid, (/(n,n=cft_lb,cft_ub)/)), subname)
       end if

       call check_ret(nf_inq_varid(ncid, 'PFTDATA_MASK', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, pftdata_mask), subname)

       call check_ret(nf_inq_varid(ncid, 'LANDFRAC_PFT', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, landfrac_pft), subname)

       call check_ret(nf_inq_varid(ncid, 'mxsoil_color', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, nsoicol), subname)

       call check_ret(nf_inq_varid(ncid, 'SOIL_COLOR', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, soicol), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_SAND', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, pctsand), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_CLAY', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, pctclay), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_WETLAND', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, pctwet), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_LAKE', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, pctlak), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_GLACIER', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, pctgla), subname)

       call check_ret(nf_inq_varid(ncid, 'GLACIER_REGION', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, glacier_region), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, pctglcmec), subname)

       call check_ret(nf_inq_varid(ncid, 'GLC_MEC', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, elevclass), subname)

       call check_ret(nf_inq_varid(ncid, 'TOPO_GLC_MEC', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, topoglcmec), subname)

       if ( outnc_3dglc )then
          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC_GIC', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglcmec_gic), subname)

          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC_ICESHEET', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglcmec_icesheet), subname)

          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_GIC', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglc_gic), subname)

          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_ICESHEET', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglc_icesheet), subname)
       end if

       if ( outnc_3dglc )then
          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC_GIC', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglcmec_gic), subname)

          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_MEC_ICESHEET', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglcmec_icesheet), subname)

          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_GIC', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglc_gic), subname)

          call check_ret(nf_inq_varid(ncid, 'PCT_GLC_ICESHEET', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, pctglc_icesheet), subname)
       end if

       call check_ret(nf_inq_varid(ncid, 'PCT_URBAN', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, urbn_classes_g), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_NATVEG', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctnatpft)), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_CROP', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft)), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctnatpft)), subname)

       if (num_cft > 0) then
          call check_ret(nf_inq_varid(ncid, 'PCT_CFT', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctcft)), subname)
       end if

       call harvdata%getFieldsIdx( harvind1D, harvind2D )
       do k = 1, harvdata%num1Dfields()
          call check_ret(nf_inq_varid(ncid, trim(mkharvest_fieldname(harvind1D(k),constant=.true.)), varid), subname)
          harvest1D => harvdata%get1DFieldPtr( harvind1D(k), output=.true. )
          call check_ret(nf_put_var_double(ncid, varid, harvest1D), subname)
       end do
       do k = 1, harvdata%num2Dfields()
          call check_ret(nf_inq_varid(ncid, trim(mkharvest_fieldname(harvind2D(k),constant=.true.)), varid), subname)
          harvest2D => harvdata%get2DFieldPtr( harvind2D(k), output=.true. )
          call check_ret(nf_put_var_double(ncid, varid, harvest2D), subname)
       end do
       deallocate( harvind1D, harvind2D )

       call check_ret(nf_inq_varid(ncid, 'FMAX', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, fmax), subname)

       call check_ret(nf_inq_varid(ncid, 'gdp', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, gdp), subname)

       call check_ret(nf_inq_varid(ncid, 'peatf', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, fpeat), subname)


       !    call check_ret(nf_inq_varid(ncid, 'Avg_Depth_Median', varid), subname)
       call check_ret(nf_inq_varid(ncid, 'zbedrock', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, soildepth), subname)

       call check_ret(nf_inq_varid(ncid, 'abm', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, agfirepkmon), subname)

       call check_ret(nf_inq_varid(ncid, 'SLOPE', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, slope), subname)

       call check_ret(nf_inq_varid(ncid, 'STD_ELEV', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, topo_stddev), subname)

       if ( outnc_vic )then
          call check_ret(nf_inq_varid(ncid, 'binfl', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, vic_binfl), subname)

          call check_ret(nf_inq_varid(ncid, 'Ws', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, vic_ws), subname)

          call check_ret(nf_inq_varid(ncid, 'Dsmax', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, vic_dsmax), subname)

          call check_ret(nf_inq_varid(ncid, 'Ds', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, vic_ds), subname)
       end if

       call check_ret(nf_inq_varid(ncid, 'LAKEDEPTH', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, lakedepth), subname)

       call check_ret(nf_inq_varid(ncid, 'EF1_BTR', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, ef1_btr), subname)

       call check_ret(nf_inq_varid(ncid, 'EF1_FET', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, ef1_fet), subname)

       call check_ret(nf_inq_varid(ncid, 'EF1_FDT', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, ef1_fdt), subname)

       call check_ret(nf_inq_varid(ncid, 'EF1_SHR', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, ef1_shr), subname)

       call check_ret(nf_inq_varid(ncid, 'EF1_GRS', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, ef1_grs), subname)

       call check_ret(nf_inq_varid(ncid, 'EF1_CRP', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, ef1_crp), subname)

       call check_ret(nf_inq_varid(ncid, 'ORGANIC', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, organic), subname)

       call check_ret(nf_inq_varid(ncid, 'URBAN_REGION_ID', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, urban_region), subname)

       ! Synchronize the disk copy of a netCDF dataset with in-memory buffers

       call check_ret(nf_sync(ncid), subname)

       ! ----------------------------------------------------------------------
       ! Make Urban Parameters from raw input data and write to surface dataset 
       ! Write to netcdf file is done inside mkurbanpar routine
       ! ----------------------------------------------------------------------

       write(6,*)'calling mkurbanpar'
       call mkurbanpar(datfname=mksrf_furban, ncido=ncid, region_o=urban_region, &
            urbn_classes_gcell_o=urbn_classes_g, &
            urban_skip_abort_on_invalid_data_check=urban_skip_abort_on_invalid_data_check)

       ! ----------------------------------------------------------------------
       ! Make LAI and SAI from 1/2 degree data and write to surface dataset 
       ! Write to netcdf file is done inside mklai routine
       ! ----------------------------------------------------------------------

       write(6,*)'calling mklai'
       call mklai(ldomain, mapfname=map_flai, datfname=mksrf_flai, &
               ndiag=ndiag, ncido=ncid )

       ! Close surface dataset

       call check_ret(nf_close(ncid), subname)

       write (6,'(72a1)') ("-",n=1,60)
       write (6,*)' land model surface data set successfully created for ', &
            'grid of size ',ns_o

    else  ! fsurdat == ' '

       write (6,*) 'fsurdat is blank: skipping writing surface dataset'

    end if  ! if (fsurdat /= ' ')

    ! Deallocate arrays NOT needed for dynamic-pft section of code

    deallocate ( organic )
    deallocate ( ef1_btr, ef1_fet, ef1_fdt, ef1_shr, ef1_grs, ef1_crp )
    deallocate ( pctglcmec, topoglcmec)
    if ( outnc_3dglc ) deallocate ( pctglc_gic, pctglc_icesheet)
    deallocate ( elevclass )
    deallocate ( fmax )
    deallocate ( pctsand, pctclay )
    deallocate ( soicol )
    deallocate ( gdp, fpeat, agfirepkmon )
    deallocate ( soildepth )
    deallocate ( topo_stddev, slope )
    deallocate ( vic_binfl, vic_ws, vic_dsmax, vic_ds )
    deallocate ( lakedepth )
    deallocate ( glacier_region )

    call harvdata%clean()

    ! ----------------------------------------------------------------------
    ! Create dynamic land use dataset if appropriate
    ! ----------------------------------------------------------------------

    if (mksrf_fdynuse /= ' ') then

       write(6,*)'creating dynamic land use dataset'

       allocate(pctlnd_pft_dyn(ns_o))
       call mkharvest_init( ns_o, spval, harvdata, mksrf_fhrvtyp )

       if (fdyndat == ' ') then
          write(6,*)' must specify fdyndat in namelist if mksrf_fdynuse is not blank'
          stop
       end if

       ! Define dimensions and global attributes

       call mkfile(ldomain, fdyndat, harvdata, dynlanduse=.true.)

       ! Write fields other pft to dynamic land use dataset

       call domain_write(ldomain, fdyndat)

       call check_ret(nf_open(trim(fdyndat), nf_write, ncid), subname)
       call check_ret(nf_set_fill (ncid, nf_nofill, omode), subname)

       call check_ret(nf_inq_varid(ncid, 'natpft', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, (/(n,n=natpft_lb,natpft_ub)/)), subname)

       if (num_cft > 0) then
          call check_ret(nf_inq_varid(ncid, 'cft', varid), subname)
          call check_ret(nf_put_var_int(ncid, varid, (/(n,n=cft_lb,cft_ub)/)), subname)
       end if

       call check_ret(nf_inq_varid(ncid, 'PFTDATA_MASK', varid), subname)
       call check_ret(nf_put_var_int(ncid, varid, pftdata_mask), subname)

       call check_ret(nf_inq_varid(ncid, 'LANDFRAC_PFT', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, landfrac_pft), subname)

       ! Synchronize the disk copy of a netCDF dataset with in-memory buffers

       call check_ret(nf_sync(ncid), subname)

       ! Read in each dynamic pft landuse dataset

       nfdyn = getavu(); call opnfil (mksrf_fdynuse, nfdyn, 'f')

       pctnatpft_max = pctnatpft
       pctcft_max = pctcft

       ntim = 0
       do 
          ! Read input pft data

          read(nfdyn, '(A195,1x,I4)', iostat=ier) string, year
          if (ier /= 0) exit
          !
          ! If pft fraction override is set, than intrepret string as PFT and harvesting override values
          !
          if ( all_veg )then
             fname = ' '
             fhrvname  = ' '
             call mkpft_parse_oride(string)
             call mkharvest_parse_oride(string)
	     write(6, '(a, i4, a)') 'PFT and harvesting values for year ', year, ' :'
             write(6, '(a, a)') '    ', trim(string)
          !
          ! Otherwise intrepret string as a filename with PFT and harvesting values in it
          !
          else
             fname = string
	     write(6,*)'input pft dynamic dataset for year ', year, ' is : ', trim(fname)
             read(nfdyn, '(A195,1x,I4)', iostat=ier) fhrvname, year2
             if ( year2 /= year ) then
                write(6,*) subname, ' error: year for harvest not equal to year for PFT files'
                call abort()
             end if
          end if
          ntim = ntim + 1

          ! Create pctpft data at model resolution
          
          call mkpft(ldomain, mapfname=map_fpft, fpft=fname, &
               ndiag=ndiag, pctlnd_o=pctlnd_pft_dyn, pctnatpft_o=pctnatpft, pctcft_o=pctcft )

          ! Create harvesting data at model resolution

          call mkharvest( ldomain, mapfname=map_fharvest, datfname=fhrvname, &
               ndiag=ndiag, harvdata=harvdata )

          ! Consistency check on input land fraction

          do n = 1,ns_o
             if (pctlnd_pft_dyn(n) /= pctlnd_pft(n)) then
                write(6,*) subname,' error: pctlnd_pft for dynamics data = ',&
                     pctlnd_pft_dyn(n), ' not equal to pctlnd_pft for surface data = ',&
                     pctlnd_pft(n),' at n= ',n
                if ( trim(fname) == ' ' )then
                   write(6,*) ' PFT string = ', string
                else
                   write(6,*) ' PFT file = ', fname
                end if
                call abort()
             end if
          end do

          call change_landuse(ldomain, dynpft=.true.)

          call normalizencheck_landuse(ldomain)
	  
          call update_max_array(pctnatpft_max,pctnatpft)
          call update_max_array(pctcft_max,pctcft)

          ! Output time-varying data for current year

          call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT', varid), subname)
          call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctnatpft))

          call check_ret(nf_inq_varid(ncid, 'PCT_CROP', varid), subname)
          call ncd_put_time_slice(ncid, varid, ntim, get_pct_l2g_array(pctcft))

          if (num_cft > 0) then
             call check_ret(nf_inq_varid(ncid, 'PCT_CFT', varid), subname)
             call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctcft))
          end if

          call harvdata%getFieldsIdx( harvind1D, harvind2D )
          do k = 1, harvdata%num1Dfields()
             call check_ret(nf_inq_varid(ncid, trim(mkharvest_fieldname(harvind1D(k),constant=.false.)), varid), subname)
             harvest1D => harvdata%get1DFieldPtr( harvind1D(k), output=.true. )
             call ncd_put_time_slice(ncid, varid, ntim, harvest1D)
          end do
          do k = 1, harvdata%num2Dfields()
             call check_ret(nf_inq_varid(ncid, trim(mkharvest_fieldname(harvind2D(k),constant=.false.)), varid), subname)
             harvest2D => harvdata%get2DFieldPtr( harvind2D(k), output=.true. )
             call ncd_put_time_slice(ncid, varid, ntim, harvest2D)
          end do
          deallocate( harvind1D, harvind2D )

          call check_ret(nf_inq_varid(ncid, 'YEAR', varid), subname)
          call check_ret(nf_put_vara_int(ncid, varid, ntim, 1, year), subname)

          call check_ret(nf_inq_varid(ncid, 'time', varid), subname)
          call check_ret(nf_put_vara_int(ncid, varid, ntim, 1, year), subname)

          call check_ret(nf_inq_varid(ncid, 'input_pftdata_filename', varid), subname)
          call check_ret(nf_put_vara_text(ncid, varid, (/ 1, ntim /), (/ len_trim(string), 1 /), trim(string) ), subname)

	  ! Synchronize the disk copy of a netCDF dataset with in-memory buffers

	  call check_ret(nf_sync(ncid), subname)

       end do   ! end of read loop

       call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT_MAX', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctnatpft_max)), subname)

       call check_ret(nf_inq_varid(ncid, 'PCT_CROP_MAX', varid), subname)
       call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)

       if (num_cft > 0) then
          call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)
          call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctcft_max)), subname)
       end if

       call check_ret(nf_close(ncid), subname)

    end if   ! end of if-create dynamic landust dataset   

    ! ----------------------------------------------------------------------
    ! Close diagnostic dataset
    ! ----------------------------------------------------------------------

    close (ndiag)
    write (6,*)
    write (6,*) 'Surface data output file = ',trim(fsurdat)
    write (6,*) '   This file contains the land model surface data'
    write (6,*) 'Diagnostic log file      = ',trim(fsurlog)
    write (6,*) '   See this file for a summary of the dataset'
    write (6,*)

    write (6,*) 'Successfully created surface dataset'

!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: change_landuse
!
! !INTERFACE:
subroutine change_landuse( ldomain, dynpft )
!
! !DESCRIPTION:
!
! Do landuse changes such as for the poles, etc.
!
! !USES:
    implicit none
!
! !ARGUMENTS:
    type(domain_type)   :: ldomain
    logical, intent(in)  :: dynpft   ! if part of the dynpft section of code

!
! !REVISION HISTORY:
! 9/10/09: Erik Kluzek spin off subroutine from original embedded code
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: n,ns_o                       ! indices
    character(len=32) :: subname = 'change_landuse'  ! subroutine name
!-----------------------------------------------------------------------

    ns_o = ldomain%ns
    do n = 1,ns_o

       ! If have pole points on grid - set south pole to glacier
       ! north pole is assumed as non-land
       
       if (abs((ldomain%latc(n) - 90._r8)) < 1.e-6_r8) then
          pctlak(n)   = 0._r8
          pctwet(n)   = 0._r8
          pcturb(n)   = 0._r8
          pctgla(n)   = 100._r8
          call pctnatpft(n)%set_pct_l2g(0._r8)
          call pctcft(n)%set_pct_l2g(0._r8)
          if ( .not. dynpft )then
             organic(n,:)   = 0._r8
             ef1_btr(n)     = 0._r8
             ef1_fet(n)     = 0._r8
             ef1_fdt(n)     = 0._r8
             ef1_shr(n)     = 0._r8
             ef1_grs(n)     = 0._r8
             ef1_crp(n)     = 0._r8
          end if
       end if

    end do

end subroutine change_landuse

!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: normalizencheck_landuse
!
! !INTERFACE:
subroutine normalizencheck_landuse(ldomain)
!
! !DESCRIPTION:
!
! Normalize land use and make sure things add up to 100% as well as
! checking that things are as they should be.
!
! Precondition: pctlak + pctwet + pcturb + pctgla <= 100 (within roundoff)
!
! !USES:
    use mkpftConstantsMod , only : baregroundindex
    use mkpftUtilsMod     , only : adjust_total_veg_area
    implicit none
! !ARGUMENTS:
    type(domain_type)   :: ldomain
!
! !REVISION HISTORY:
! 9/10/09: Erik Kluzek spin off subroutine from original embedded code
!
!EOP
!
! !LOCAL VARIABLES:
    integer  :: m,k,n,ns_o                  ! indices
    integer  :: nsmall                      ! number of small PFT values for a single check
    integer  :: nsmall_tot                  ! total number of small PFT values in all grid cells
    real(r8) :: suma                        ! sum for error check
    real(r8) :: suma2                       ! another sum for error check
    real(r8) :: new_total_veg_pct           ! new % veg (% of grid cell, total of natural veg & crop)
    real(r8) :: bare_pct_p2g                ! % of bare soil, as % of grid cell
    real(r8) :: bare_urb_diff               ! difference between bare soil and urban %
    real(r8) :: pcturb_excess               ! excess urban % not accounted for by bare soil
    real(r8) :: sum8, sum8a                 ! sum for error check
    real(r4) :: sum4a                       ! sum for error check
    real(r8), parameter :: tol_loose = 1.e-4_r8               ! tolerance for some 'loose' error checks
    real(r8), parameter :: toosmallPFT = 1.e-10_r8            ! tolerance for PFT's to ignore
    character(len=32) :: subname = 'normalizencheck_landuse'  ! subroutine name
!-----------------------------------------------------------------------

    ! ------------------------------------------------------------------------
    ! Normalize vegetated area so that vegetated + special area is 100%
    ! ------------------------------------------------------------------------

    ns_o = ldomain%ns
    do n = 1,ns_o

       ! Check preconditions
       if ( pctlak(n) < 0.0_r8 )then
          write(6,*) subname, ' ERROR: pctlak is negative!'
          write(6,*) 'n, pctlak = ', n, pctlak(n)
          call abort()
       end if
       if ( pctwet(n) < 0.0_r8 )then
          write(6,*) subname, ' ERROR: pctwet is negative!'
          write(6,*) 'n, pctwet = ', n, pctwet(n)
          call abort()
       end if
       if ( pcturb(n) < 0.0_r8 )then
          write(6,*) subname, ' ERROR: pcturb is negative!'
          write(6,*) 'n, pcturb = ', n, pcturb(n)
          call abort()
       end if
       if ( pctgla(n) < 0.0_r8 )then
          write(6,*) subname, ' ERROR: pctgla is negative!'
          write(6,*) 'n, pctgla = ', n, pctgla(n)
          call abort()
       end if

       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
       if (suma > (100._r8 + tol_loose)) then
          write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla must be'
          write(6,*) '<= 100% before calling this subroutine'
          write(6,*) 'n, pctlak, pctwet, pcturb, pctgla = ', &
               n, pctlak(n), pctwet(n), pcturb(n), pctgla(n)
          call abort()
       end if

       ! First normalize vegetated (natural veg + crop) cover so that the total of
       ! (vegetated + (special excluding urban)) is 100%. We'll deal with urban later.
       !
       ! Note that, in practice, the total area of natural veg + crop is typically 100%
       ! going into this routine. However, the following code does NOT rely on this, and
       ! will work properly regardless of the initial area of natural veg + crop (even if
       ! that initial area is 0%).
       
       suma = pctlak(n)+pctwet(n)+pctgla(n)
       new_total_veg_pct = 100._r8 - suma
       ! correct for rounding error:
       new_total_veg_pct = max(new_total_veg_pct, 0._r8)

       call adjust_total_veg_area(new_total_veg_pct, pctnatpft=pctnatpft(n), pctcft=pctcft(n))

       ! Make sure we did the above rescaling correctly

       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
       if (abs(suma - 100._r8) > tol_loose) then
          write(6,*) subname, ' ERROR in rescaling veg based on (special excluding urban'
          write(6,*) 'suma = ', suma
          call abort()
       end if

       ! Now decrease the vegetated area to account for urban area. Urban needs to be
       ! handled specially because we replace bare soil preferentially with urban, rather
       ! than rescaling all PFTs equally.

       if (pcturb(n) > 0._r8) then
          
          ! Replace bare soil preferentially with urban
          bare_pct_p2g = pctnatpft(n)%get_one_pct_p2g(baregroundindex)
          bare_urb_diff = bare_pct_p2g - pcturb(n)
          bare_pct_p2g = max(0._r8, bare_urb_diff)
          call pctnatpft(n)%set_one_pct_p2g(baregroundindex, bare_pct_p2g)
          pcturb_excess = abs(min(0._r8,bare_urb_diff))
          
          ! For any urban not accounted for by bare soil, replace other PFTs
          ! proportionally
          if (pcturb_excess > 0._r8) then
             ! Note that, in this case, we will have already reduced bare ground to 0%

             new_total_veg_pct = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() - pcturb_excess
             if (new_total_veg_pct < 0._r8) then
                if (abs(new_total_veg_pct) < tol_loose) then
                   ! only slightly less than 0; correct it
                   new_total_veg_pct = 0._r8
                else
                   write(6,*) subname, ' ERROR: trying to replace veg with urban,'
                   write(6,*) 'but pcturb_excess exceeds current vegetation percent'
                   call abort()
                end if
             end if

             call adjust_total_veg_area(new_total_veg_pct, pctnatpft=pctnatpft(n), pctcft=pctcft(n))
          end if

       end if ! pcturb(n) > 0

       ! Confirm that we have done the rescaling correctly: now the sum of all landunits
       ! should be 100%
       suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n)
       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
       if (abs(suma - 100._r8) > tol_loose) then
          write(6,*) subname, ' ERROR: landunits do not sum to 100%'
          write(6,*) 'n, suma, pctlak, pctwet, pctgla, pcturb, pctnatveg, pctcrop = '
          write(6,*) n, suma, pctlak(n), pctwet(n), pctgla(n), pcturb(n), &
               pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g()
          call abort()
       end if
          
    end do

    ! ------------------------------------------------------------------------
    ! Do other corrections and error checks
    ! ------------------------------------------------------------------------

    nsmall_tot = 0

    do n = 1,ns_o
       
       ! If the coverage of any PFT or CFT is too small at the gridcell level, set its
       ! % cover to 0, then renormalize everything else as needed
       call pctnatpft(n)%remove_small_cover(toosmallPFT, nsmall)
       nsmall_tot = nsmall_tot + nsmall
       call pctcft(n)%remove_small_cover(toosmallPFT, nsmall)
       nsmall_tot = nsmall_tot + nsmall

       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
       if ( abs(suma - 100.0_r8) > 2.0*epsilon(suma) )then
          pctlak(n)    = pctlak(n)    * 100._r8/suma
          pctwet(n)    = pctwet(n)    * 100._r8/suma
          pcturb(n)    = pcturb(n)    * 100._r8/suma
          pctgla(n)    = pctgla(n)    * 100._r8/suma
          call pctnatpft(n)%set_pct_l2g(pctnatpft(n)%get_pct_l2g() * 100._r8/suma)
          call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma)
       end if
       
       ! Roundoff error fix
       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
       suma2 = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
       if ( (suma < 100._r8 .and. suma > (100._r8 - 1.e-6_r8)) .or. &
            (suma2 > 0.0_r8 .and. suma2 <  1.e-6_r8) ) then
          write (6,*) 'Special land units near 100%, but not quite for n,suma =',n,suma
          write (6,*) 'Adjusting special land units to 100%'
          if (pctlak(n) >= 25._r8) then
             pctlak(n) = 100._r8 - (pctwet(n) + pcturb(n) + pctgla(n))
          else if (pctwet(n) >= 25._r8) then
             pctwet(n) = 100._r8 - (pctlak(n) + pcturb(n) + pctgla(n))
          else if (pcturb(n) >= 25._r8) then
             pcturb(n) = 100._r8 - (pctlak(n) + pctwet(n) + pctgla(n))
          else if (pctgla(n) >= 25._r8) then
             pctgla(n) = 100._r8 - (pctlak(n) + pctwet(n) + pcturb(n))
          else
             write (6,*) subname, 'Error: sum of special land units nearly 100% but none is >= 25% at ', &
                  'n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n),suma = ', &
                  n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
                  pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(),suma
             call abort()
          end if
          call pctnatpft(n)%set_pct_l2g(0._r8)
          call pctcft(n)%set_pct_l2g(0._r8)
       end if
       if ( any(pctnatpft(n)%get_pct_p2g() > 0.0_r8 .and. pctnatpft(n)%get_pct_p2g() < toosmallPFT ) .or. &
            any(pctcft(n)%get_pct_p2g()    > 0.0_r8 .and. pctcft(n)%get_pct_p2g()    < toosmallPFT )) then
          write (6,*) 'pctnatpft or pctcft is small at n=', n
          write (6,*) 'pctnatpft%pct_p2l = ', pctnatpft(n)%get_pct_p2l()
          write (6,*) 'pctcft%pct_p2l = ', pctcft(n)%get_pct_p2l()
          write (6,*) 'pctnatpft%pct_l2g = ', pctnatpft(n)%get_pct_l2g()
          write (6,*) 'pctcft%pct_l2g = ', pctcft(n)%get_pct_l2g()
          call abort()
       end if
       
       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
       if (suma < 100._r8-epsilon(suma) .and. suma > (100._r8 - 4._r8*epsilon(suma))) then
          write (6,*) subname, 'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', &
               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
               pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g()
          call abort()
       end if
       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
       if ( abs(suma-100._r8) > 1.e-10_r8) then
          write (6,*) subname, ' error: sum of pctlak, pctwet,', &
               'pcturb, pctgla, pctnatveg and pctcrop is NOT equal to 100'
          write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,sum= ', &
               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
               pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(), suma
          call abort()
       end if
       
    end do

    ! Check that when pctnatveg+pctcrop identically zero, sum of special landunits is identically 100%

    if ( .not. outnc_double )then
       do n = 1,ns_o
          sum8  =         real(pctlak(n),r4)
          sum8  = sum8  + real(pctwet(n),r4)
          sum8  = sum8  + real(pcturb(n),r4)
          sum8  = sum8  + real(pctgla(n),r4)
          sum4a =         real(pctnatpft(n)%get_pct_l2g(),r4)
          sum4a = sum4a + real(pctcft(n)%get_pct_l2g(),r4)
          if ( sum4a==0.0_r4 .and. sum8 < 100._r4-2._r4*epsilon(sum4a) )then
             write (6,*) subname, ' error: sum of pctlak, pctwet,', &
                  'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8
             write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', &
                  n,pctlak(n),pctwet(n),pcturb(n),pctgla(n), &
                  pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g()
             call abort()
          end if
       end do
    else
       do n = 1,ns_o
          sum8  =         pctlak(n)
          sum8  = sum8  + pctwet(n)
          sum8  = sum8  + pcturb(n)
          sum8  = sum8  + pctgla(n)
          sum8a =         pctnatpft(n)%get_pct_l2g()
          sum8a = sum8a + pctcft(n)%get_pct_l2g()
          if ( sum8a==0._r8 .and. sum8 < (100._r8-4._r8*epsilon(sum8)) )then
             write (6,*) subname, ' error: sum of pctlak, pctwet,', &
                  'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8
             write (6,*) 'Total error, error/epsilon = ',100._r8-sum8, ((100._r8-sum8)/epsilon(sum8))
             write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,epsilon= ', &
                  n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
                  pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(), epsilon(sum8)
             call abort()
          end if
       end do
    end if

    ! Make sure that there is no vegetation outside the pft mask
    do n = 1,ns_o
       if (pftdata_mask(n) == 0 .and. (pctnatpft(n)%get_pct_l2g() > 0 .or. pctcft(n)%get_pct_l2g() > 0)) then
          write (6,*)'vegetation found outside the pft mask at n=',n
          write (6,*)'pctnatveg,pctcrop=', pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g()
          call abort()
       end if
    end do

    ! Make sure that sums at the landunit level all add to 100%
    ! (Note that we don't check pctglcmec here, because it isn't computed at the point
    ! that this subroutine is called -- but the check of sum(pctglcmec) is done in
    ! mkglcmecMod)
    ! (Also note that we don't need to check pctnatpft or pctcft, because a similar check
    ! is done internally by the pct_pft_type routines.)
    do n = 1,ns_o
       if (abs(sum(urbn_classes(n,:)) - 100._r8) > 1.e-12_r8) then
          write(6,*) 'sum(urbn_classes(n,:)) != 100: ', n, sum(urbn_classes(n,:))
          call abort()
       end if
    end do

    if ( nsmall_tot > 0 )then
       write (6,*)'number of small pft = ', nsmall_tot
    end if

end subroutine normalizencheck_landuse

end program mksurfdat
